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Summary 


The application of the numerical method of characteristics to the solution 
of the differential equations defining unsteady flow in partially filled 
drainage system sized pipes is outlined. 


The derivation of the flow equations is presented, together with the neces- 
Sary boundary equation formulation to represent variable inflow, system 
discharge and leakage flow past a stationary deposited solid. 


A computer program, written in Fortran, is included, together with typical 
output, that establishes the applicability of this computational method tu 
unsteady flow analysis in gravity flow drainage systems. 


Proposals for the extension of the described techniques to the prediction 
of solid transport and flow attenuation in long pipes are also presented. 
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This report is one of a group documenting National Bureau of Standards 
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Urban Development /Uffice of Policy Development and Research, Division of 
Energy Building Technology and Standards, under HUD Interagency Agreement 
H-48-78. 


Report prepared by Dr. Je A. Swaftield, Senior Lecturer, Drainage Research 
Group, Department of Building Technology, Brunel University, Uxbridge, Uk., 
during a study leave period as a guest research worker at NBS/Stevens 
Institute of Technology. 
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Experimenta] results included in this report are drawn from the published wor 
of the Drainage Research Group at Brunel University. 
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l. INTRODUCTION 


Unsteady flow is defined as fluid flow wherein the flow parameters vary with 
time and as such will occur in all fluid carrying conduits if the system 
boundary conditions are rapidly changed, either by design or as a result of 
some unforeseen event. The rate of change of the system boundary conditions 
determine the severity of the flow parameter changes that result. In many 
cases the resulting pressures generated in the flow system determine the 
design conditions. 


The most well documented area of transient analysis refers to the propagation 
of pressure surges through full bore flow carrying pipe systems. For example 
the pressure variations following inadvertent power failure to a pipeline 
pumping station often provides the system design condition and leads to the 
need to introduce safety systems such as air chambers, increased pump inertia 
or by pass systems. Similarly the surges generated by load rejection by hydro- 
electric power station turbines lead to surge propagation in the water supply 
tunnels, normally dealt with by the inclusion of surge shafts cut into the 

rock surrounding the hydroplant supply tunnels. 


The vast majority of research into unsteady flow has therefore been directed 
towards the large civil engineering applications such as power stations, cross 
country oil, water or sewage pipelines. Currently these areas are well docu- 
mented and several texts are available describing appropriate analysis and 
design techniques [1,2]. 


The prediction techniques currently employed may be traced back to the late 
19th century and for many years were based on the Bergeron graphical methods, 
that incidentally also find application in electrical surge problems. 


The widespread availability of digital computers from the 1960's introduced or 
made practical, more versatile methods. Although translations of the graphi- 
cal techniques were used initially, the numerical method of characteristics 
employed to solve the differential equations defining the unsteady flow has 
now come to be widely accepted as the more versatile technique [3]. 


As mentioned the applications to be found in the literature refer almost 
exclusively to large scale civil engineering applications. However, the prob- 
lem of unsteady flow translates across the boundaries of scale, being dependent 
on the relationship between the rate of change of boundary conditions and the 
propagation rate of transient pressures within the system. Such applications 
may be found in a wide range of small scale systems, including, for example, 
diesel fuel injection systems. The analysis technique based on the method of 
characteristics solution was successfully applied to a transient analysis of 
the Concorde fuel system, both in the ground refuelling mode and the in flight 
fuel transfer and emergency ejection modes [4]. Similarly the methodology has 
been applied to airport hydrant refuelling networks. 


Although this vast bank of experience and literature exists it is, as 
mentioned, primarily directed to full bore flow large civil engineering 
applications. 


Unsteady flow may, by definition, also occur in partially filled pipe or open 
channel flow. Again the available literature applies to the large excavated 
straight sided channels to be found in such applications as power station tur- 
bine discharge channels. In the open channel flow situation the full bore flow 
transient pressure changes are replaced by channel depth variations, pressure 
waves moving at the appropriate sonic velocity being replaced by surface waves. 
The simplification possible in air or vapor free fluid full bore flow that the 
wave propagation speed is constant, depending only on the liquid and pipe 
materials and dimensions, no longer holds as the wave propagation speed depends 
on flow depth and channel shape. Thus an extra variable is effectively added 
to the local flow velocity and pressure, or depth, namely the local wave speed. 


As will be seen in the analysis presented this also introduces other difficul- 
ties in the numerical solution of the governing equations since the uniform x-t 
grid system generally employed has to be modified to allow for variable charac- 
teristic slope as the local wave speed varies. 


As with the scale translation from hydroelectric plant to aircraft fuel system, 
so too does the basic flow analysis translate from large excavated open chan- 
nels to partially filled drain piping systems. The boundary conditions change 
to represent, for example, varying inflow dependent on the appliances connected 
to the system. The channel cross sectional shape does not constrain the appli- 
cation of the fundamental analytical approach. 


At present no numerical analysis technique has been documented to predict the 
depth and flow rate along gravity driven unsteady flow in partially filled 
drainage sized systems; however the analysis techniques incorporating the 
method of characteristics together with boundary conditions expressed in terms 
of time dependent depth or flow rates may be employed. 


The current report is designed to act as an introduction to the use of these 
methods to predict depth-time profiles along a simulated drain. In particular 
the analysis will attempt to predict the depth history upstream of a deposited 
solid in the drain pipe as a prelude to a consideration of the forces acting on 
the solid and its subsequent motion along the drain pipe. 


In addition the analysis techniques are shown to be applicable to determine 
flow attenuation in long drain pipes, this topic being considered increasingly 
important due to the probable reduction in overall system flow rates due to 
water conservation proposals. 


Similarly the techniques can be applied to pipe sizing calculations in order to 
avoid the occurrence of local full bore flow and associated venting problems. 


The differential equations defining unsteady flow are presented in this report, 
together with the application of the method of characteristics to yield a total 
differential equation which may be solved by numerical methods. The method of 
specific time intervals is presented together with the necessary formulation of 
boundary equations to represent variable inflow, pipe discharge and leakage 
flow past the solid. Finally a computer program designed to analyze the trans- 
ient depth response to a stationary solid in a pipe supplied with a variable 
inflow is presented together with typical output for flow and system variables 
representative of drainage design values. ; 


Extension of the methods described to flow attenuation and solid motion 
prediction are also proposed. 

Le UNSTEADY FLOW IN OPEN CHANNELS 

2.1 DEFINING EQUATIONS 

The equations defining one dimensional unsteady flow are presented for open 
channels. Frictional losses are assumed to be adequately represented by the 
Manning equation utilizing the local velocity and flow depth. Channel slopes 


are assumed small enough for cosine slope (cosa) to be approximately 1. 


Figure 1 illustrates an elemental fluid strip. Application of the unsteady 
momentum equation yields: 


(i) Net hydrostatic force: pg dh Ax A 
dx 
(ii) Shear force on wetted: ese 


perimeter P 


(iii) Gravity force in flow: pg A Ax sina 
direction 
(iv) Net efflux momentum: 2 (pv2A) Ax 
ox 
(v) Increase of momentum within: 0 (pAV) Ax 
ot 


Thus the momentum equation becomes 


—pgAxA a - TP Ax + pgAAx sina 
x 
(1) 
= 9 (pv2A) Ax + 2_(pAV) Ax 
ox ot 


Expanding and dividing by pAax 
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DEIN ris pried port ghia str lop let 
pees nate) sina wo meee 
ax om 2 ax | A 3x 
2) 
+ Vga fav. oO 
Avot 28 & 
where m = A/P 
The continuity equation applied to the control volume of Figure 1 yields 
—a(pAV) AX = 3 (pAax) (3) 
Ox at 
Expanding and dividing by pAx 
vy 9A +2444 2V 20 (4) 
aX at ax 
Multiplying (4) by V/A and subtracting from equation (2) yields 
g2h + to - gsim + ¥V 4 aV = 9 ©) 
aX om aX t 
Let So = sina and gS = To/pm 
where S is the slope of the energy grade line as defined by the Manning 
equation 
es stele (0) 
nA] 3 
Note that S = S, only under steady uniform flow conditions. 
Equation 5 becomes 
ah + eS — S,) + VON + ot =o (7) 
esac ar Oo iran) ae 


From Figure 1 9A = dh. T. 


where T is surface width. 


Hence — = psi 
ax ax 
dA dh 

nd ie) tea 

3 at at 


ae Leong Cid oe OM Ag (8) 
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Thus equations 7 and 8 represent the unsteady flow conditions in terms of 
local depth and average velocity. These equations may be solved by means of 
the method of characteristics. 

2.2 SOLUTION BY METHOD OF CHARACTERISTICS 

Equation 7 and 8 may be combined as 

where L, is equation 7 and Ly equation 8 


Oy eA ee aor ea (V +8) + Sh 
ax at ox d t 


] 2 (2) 
ap pis) Sad = © 
In order to solve equations 7 and 8 it is necessary to transform the equations 
into a total derivative expression. 


For the terms in bracket [1] to become the total derivative dV/dt, it is 
necessary for 


V +A = 9X 
dt 


and for the terms in bracket [2] to become the total derivative dh/dt then 


hence V + AA = V + g/dT 


therefore A = # Re 


Qa./Q 
(anal fre 
I+ 
» 


and Ase Ea 


The term VEE has dimensions of velocity and is identified as the local wave 
AG 


speed 


Suen 
c = (10) 


aVitedh+y(s-s) =0 (11) 
dt c dt fe) 
dx c (12) 


subject toc—— =. V = 
dt 


Referring to Figure 2, if the variables V and h are known at R and S then four 
equations may be written in terms of the unknowns at point P, namely 


Bee t 
Vel ee f + dh + if g(S-S,)dt = 0 (13) 
(es 
cS hR tk 
C 
Cp 
Xp 7 Xp = iT (V+c )dt (14) 
CR 
hp 1 tp ’ 
Vp - Vote = dh + 5 g(S-S,)dt = 0 (15) 
Rao te 
- S S 
C 
Cp 
Xp =F Xo = (V-c )dt (16) 


It is stressed that these equations are paired, i.e. equation 13 only holds if 
equation 14 is satisfied and 15 only applied if 16 is satisfied. Generally 
these are referred to as ct and C~ characteristics as shown above. 


In most cases~a first order integration is satisfactory, however attention must 
be paid to the choice of time step to ensure a stable solution, as mentioned by 
Fox [1]. 


In terms of Figure 2 the time step must be sufficiently small to ensure that 
points R and S fall within + Ax on either side of point P. 


Generally a suitable rule is to set 

At = Ax/2c (17) 
where c is the wave speed appropriate to the initial flow in the channel prior 
to the initiation of unsteady flow. As will be seen the solution requires an 


initial steady flow in the channel, although this may be infinitesimal. 


Applying a first order approximation to equations 13 to 16, in terms of 
Figure 2 yields 


pane pot a (hei= tet e (Se St RCS = 0 (18) 
Xp — Xp = (Vp + cp) At (19) 
VE Ve = Gis note (Sy — Sati = 0 (20) 
Kp) = xce— (Wig =e) Mt : (21) 


It will be noted that conditions at R and S at time (t-At) are determined by 
interpolation between points A and C and C and B respectively. This interpola- 
tion introduces a dispersive element to the calculation as effects arriving at 
A, C and B at time t - At are assumed, via the interpolation, to effect con- 
ditions at R and S at time t - At. This is comparable to increasing the 

local wave speed. 


Interpolation effects should be minimized by arranging for R and S to fall as 
close as possible to A and B by adjusting the time step At. 


Two flow regimes may be identified for open channels defined in terms of the 
Froude Number Fr = V/vgh 


1)  Subcritical flow, Froude N° < 1 


Here the local wave speed is greater than the flow average velocity. 
Thus waves may be propagated upstream. 


2) Supercritical flow, Froude N° > 1 


Here the local wave speed is less than the average flow velocity at 
that section, hence waves may not be propagated upstream. Flow may 
be transformed from supercritical to subcritical via a hydraulic 
jump. 


Figure 3 illustrates the importance of these two flow regimes on the solution 
of equations 13 to 16. 


If c > V then the conditions at P are determined by the intersection of the 
ct and C7 drawn from P into the AC and BC sections. 


If c < V then conditions in the downstream section BC cannot effect point P. 
The slope of the C” characteristic, PS, becomes positive and both R and § lie 
in the AC section as shown. 


For the subcritical flow encountered in the drainage applications considered 
and the equations derived below refer to this flow regime. Similar derivations 


may be undertaken for supercritical flow. 


Referring to Figure 2 for subcritical flow: 


Vao= Vi - x 
¢ 7 Vay 4e % XA Ax 
Ca Se =X 
c ~ R = “Cc ~ R = (Ve + CR) At 
Co CA Xo Xn Ax 
ha > h 

and C R = (Vp + CR) Atl 


as Xp = Xc and Xp — Xp = (Vp + cp) At 


Solution yields 


Paper Chie Cy EY CSAN EC o AY (22) 
= R We) (CVG = Volt co Sae,) 


mice Chis VRie) )) t cavp © (2s 
R ] + Cc8 = CaO 


= Ye ee Che ad ha) (0 (Vp i CR)) (24) 


=r 
v0] 
l 


Similarly 


Ve i @(Vccp st ccVRp) 


ea oF Wen 2 e ; 
ee eminence 1B) (26) 
S hea? 6) (eS) Gin) 
he = he ate (Vs = Co) (he = hp) (27) 


The determination of conditions at P at time t + At requires the following 


steps: 


Gi) 
(ii) 


Girins) 


(iv) 


(v) 


All conditions known at time t for nodal points A BC etc. 


Values of V, h and c at interpolation points RK and S calculated fron 
equations 21 - 27. 


Using these values of V, h, and c the conditions at P, iee. velocity 
V and depth h, at time t + At are calculated by means of equations 
18 and 20. 


The value of wave speed c at P at time t + at is calculated from 
equation 10. The value of flow surface width and cross sectional 
area are calculated from flow depth, h, and the channel shape 
relationships. 


The sequence is repeated at each time step. 


2.3 APPLICATION OF SOLUTION TO DRAINAGE FLOW 


Equations 18 to 21 may be expressed as: 


C (28) 
Xpeeee Xpe (VR + CpdAt 


C AD) 
Xp = Xs = (Vo =. Co )At 


where “XI'= tgi/icp 


X3 = g/cc 

X2 = Vp + 9 hp - g(S, = Sat 
eR 

Kae a Rega s anes ete g(S< = Sat 
es 


/ 


Figure 4 illustrates a typical drainage pipe length to be analyzed in terms 
of flow depth and velocity at each section. 


The application may be dealt with in two distinct sections: 
(i) Internal or nodal points. 
The values of h, V and c at all points Ax apart between x = Ax and 
x = (L - 4x) may be calculated by the sequence set out above, namely 
by simultaneous solution of equations 28, 29 and use of equation 1V to 
yield wave speed. 


(ii) Boundary conditions. 


In order to predict h, V and c at the system boundary it is necessary 
to solve either 28 or 29 with an appropriate boundary condition. 


At pipe entry a suitable boundary condition would be the inflow profile 
Q = f(r) 
to be solved with the C. characteristic: 
Qe) = ViAy sve) 
V,; = X4 + X3 hy 
Q(t) = An (x4 X30) 
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where Aj = f (hy ) 

Q(t) = £(h,)(X4 + X3 hy) 
In the form  . 

O= 0b) = £Ghp)) €X4: =) X3%h})) (30) 
this relationship may be solved by bisection method as values of Q are known 
at each time step. The channel cross section shape relationship is also 
required so that no direct solution of equation 30 is available for circular 


section channels. 


At pipe exit in the subcritical flow regime the flow depth approaches the 
critical depth value, given by zero value of the expression: 


et Terit = 1, 
8 Acrit 

where A and T are functions of depth, h. 

This condition may be solved with the ct characteristic 
WigRE) 9 S20 a 

where N = N° of pipe length sections, each of length Ax. 


The boundary condition becomes 


Die Sigman? 
[(X2 - Xl bys) Angy ]° Tye /8 ae 1=0 (31) 


Solution may again be achieved by use of the bisection method together with the 
use of the area to depth relationship for the channel. 


2.4 APPLICATION OT WASTE SOLID BOUNDARY CONDITION 


Considering a stationary solid deposited at some point along the waste pipe, 
the water depth and velocity upstream of the solid may be predicted if a suit- 
able boundary equation may be written linking flow past the solid to upstream 
conditions. 


Figure 5 illustrates the relationship between flow past a stationary solid 

and the specific energy upstream. These results were compiled during a Brunel 
University Drainage Research Group study of solid transport in drainage 
systems. 
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From the experimental work summarized by Figure 5, the flow past the solid may 
be expressed by the following relationship: 


2 2 
Q_= K(h + one SE,) (32) 


where SE = h + v2/2, flow specific energy and SE, is the flow specific energy 
required for flow initiation past the solid (when- SEo # 0) 


Equation 32 may then be solved with the C* characteristic 
Vy+1 = X2 - Xl hy+1 
where Q = Vy+) An+] 


so that 
1 a 2 2 
Ay+, (X2 - X) byyy) = K [hyyy Sor (X2 - Xl hyy))~ 7 SEQ] 
._ This expression results in a quartic in terms of water depth upstream of the 


solid, hy+,, see Appendix 2. 


This quartic must be solved by an iterative technique as the flow area, Ay+, 
is a function of hy4). 


The Newton-Raphson method may be employed to carry out the necessary iteration 
solution. 


Once the value of hy+, has been determined the value of Vy+, and cy+, may be 
determined from equation 28 and 10. 


As mentioned the SE, term is the flow specific energy required to initiate 

flow past the stationary solid. If the value of flow specific energy at tine 
t is less that SE, then the value of flow velocity at the solid at time t + At 
is set equal to zero. The flow depth then comes directly from equation 28 as: 


hy+1 = X2/X1 (33) 


This implies that the flow depth upstream of the solid must rise to SE, prior 
to the tnitiation of flow past the solid. 


This solution is set out in detail in Appendtx 2. 


In this model no account is taken of flow downstream of the solid. 
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From Figure 5 the flow past the solid may alternatively be expressed as 


where SEy SEp are the specific energy values immediately upstream and down- 
stream of the solid. 


This boundary equation may be used to link the flow conditions upstream and 
downstream of the solid, for example the upstream flow is governed by the er 
equation 28 and that downstream is governed by the C equation 29. 


This formulation of the boundary condition would allow a series of nodal 
points downstream of the solid to be dealt with, the pipe being then term- 
inated by the boundary condition already described, i-e. critical depth at 
discharge. 


The techniques described above have been included in a Fortran program 
TRANSCA run on the NBS Center for Building Technology Perkin Elmer 732 
computer. 


This program is included in Appendix 1 together with a flow chart and 
description and sample input data while representative program output are 
included in Appendix 3. 


3. PROGRAM OUPUT AND DEVELOPMENT 


Figures 6, 7 and 8 illustrate the predicted depth and flowrate at sections 
along the simulated drain pipe considered, while Appendix 3 presents typical 
program output. 


The program calculation technique requires a known steady uniform flow to be 

set up along the channel prior to the calculation of unsteady conditions. 

In the model presented this is assumed to be established with no solid in the 
pipe- The calculation then assumes the presence of the solid from the first 

time step onwards. This explains the depth increasing wave that is shown in 

Figure 8 moving upstream from the solid boundary. This is merely a function 

of the initial conditions chosen. 


An alternative base condition would be either steady flow with the solid in 
place, with the associated back water profile upstream of the stationary 
solid, or stationary fluid trapped behind the solid. In both cases the fluid 


depth would increase at the solid location, forming the expected backwater 
curve profile. 


The former set of flow conditions could have been achieved in the current 
simulation by holding the inflow rate constant, however this would have 
required a longer computer run. The establishment of steady conditions is 
demonstrated in Figure 6 as the inflow profile is assumed to become constant. 
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The results demonstrate that the techniques described may be employed to 

yield values of depth and flowrate under unsteady flow conditions in open 
channel flow. In particular the results indicate that the concept and use of 
specific energy to leakage flow characteristic, rather akin to a valve 
pressure flow relationship, is capable of accurately describing the boundary 
conditions at a stationary solid. The values of depth at the solid with 

time obviously require experimental verification, however the depths predicted 
are generally in line with observations of water build up behind stationary 
solids made during the Brunel University drainage investigations. 


Two main development paths may be proposed based on the analysis techniques 
described: 


(i) Application to the moving solid case. 


This case would require the introduction of flow relative velocity 
into the boundary conditions describing the solid described above. It 
is proposed that, once this has been done, the stationary solid speci- 
fic energy differential across the blockage to flow past relationship 
is assumed to apply to the moving case. 


A number of unknowns are required to obtain a solution, namely the 
force balance necessary to be overcome to initiate motion and the 
solid-flow momentum equation that would determine solid velocity and 
distance traversed in each time step. 


One simplifying assumption that could be made would be to assume that 
the flow ahead of the solid, i-e. downstream, assumes the normal depth 
and velocity characteristic to the channel dimensions and the flow 
past the solid. 


This case is currently being considered. 


(ii) Application of the techniques described to the prediction of flow 
attenuation in long channels. 


At first sight it would appear that the existing program would be 
capable of dealing with this case, however the errors due to inter- 
polation made necessary by the adoption of a fixed time-distance 
grid, Figure 2, could introduce errors. The interpolation techniques 
effectively disperse the moving wave fronts and would tend to over- 
estimate the attenuation produced. 


A free grid in which both time and distance become calculated 
variables is a solution to this problem. This is not normally done as 
it results in a rather “untidy” output, however this method will be 
investigated as part of the current investigation. 
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As mentioned the method of characteristics requires that there be an initial 
steady flow condition along the channel, or that the flow is at rest at 

known depth. This limiting condition is necessary to provide the base values 
of depth, velocity and wave speed at each calculating position. The influ- 
ence of magnitude of this initial assumed flow on the attenuation of a con- 
stant superimposed inflow will also be considered as part of (ii) above, 
together with an assessment of the effect of the interpolation ratio employed 
in choosing the appropriate calculation time step. 


4. CONCLUSIONS 


The program presented indicates that the techniques involved in the solution 
of unsteady flow problems via the method of characteristics may be applied 
to partially filled drainage pipeflow. 


The methods presented will be employed to investigate further their applica- 
tion to solid waste transport analysis and to the study of flow attenuation 
in long channels, typically representative of underground building drainage 
systems. 
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Figure 8. Depth proriles along the 5m simulated drain at 14 time intervals 
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Appendix i, : 


Program TRANSCA 


This appendix presents a complete print out of this program, written in 
Fortran, together with a flow chart and sample input data. The program was 
run on the NBS Center for Building Technology Perkin Elmer 732 computer. 


It was found necessary to employ double precision calculating techniques in 
the program due to the inclusion of calculations involving the root of small 
differences between large numbers. Otherwise no special facilities are 
required. 


The program -ccepts data in SI units with the exception of the inflow 
profile, this is read in as litres/second and corrected to m/s within the 
program. 
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